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We report structural evidence of dynamic reorganization in vortex matter in clean NbSe 2 by joint 
small angle neutron scattering and ac-susceptibility measurements. The application of oscillatory 
forces in a transitional region near the order-disorder transition results in robust bulk vortex lattice 
configurations with an intermediate degree of disorder. These dynamically-originated configurations 
correlate with intermediate pinning responses previously observed, resolving a long standing debate 
regarding the origin of such responses. 

PACS numbers: 74.25.Uv, 61.05.fg, 64.60.Cn 


In a wide variety of complex systems, competing in¬ 
teractions promote an order-disorder transition (ODT). 
Ordered phases are characterized by spatial correlations 
decaying weakly over distances larger than the relevant 
system scale, whereas disordered configurations are char¬ 
acterized by correlation lengths ( of the order of the mean 
inter-particle distance ao. Configurations with intermedi¬ 
ate degrees of disorder, strong enough to affect the system 
response, but still with £ ao, have received attention 
recently [I]. Vortex matter in superconductors provides 
an ideal model system for the experimental study of the 
topic PJ: vortex-vortex interactions favoring an ordered 
vortex lattice (VL) compete with both thermal fluctua¬ 
tions and pinning interactions that tend to disorder the 
system. 

In very low-pinning superconductors, such as clean 
NbSe 2 single crystals, most of the vortex field-temper¬ 
ature phase diagram is properly described by an ordered 
dislocation-free Bragg Glass (BG) phase [3H5], which un¬ 
dergoes an ODT near the normal-superconductor transi¬ 
tion In practice, when a superconductor is cooled 

from the normal state in an external magnetic field (field- 
cooled /FC), energy barriers may trap the VL in highly 
disordered metastable configurations [8j. Even so, high 
transport current densities 0 [10] or large oscillatory 
“shaking” magnetic fields mm may anneal the VL 
into the ordered low-tenrperature BG. Once in the BG, 
as temperature or magnetic field is increased the VL soft¬ 
ens and accommodates to the random pinning potential 
more efficiently. Eventually, vortex entanglement and the 
proliferation of VL dislocations increase the effective pin¬ 
ning, producing a sudden rise of the critical current J c 
known as the Peak Effect (PE) [ 7 } [ 5 ] , which is the finger¬ 
print of the ODT in vortex matter. While both the high- 
pinning disordered phase and the low-pinning ordered 
phase, above and below the ODT, are widely accepted in 
the literature, intermediate responses have been ascribed 


to surface contamination and partial reordering by the 
probing transport current mmmm- However, in the 
narrow transitional region adjacent to the PE, for which a 
multidomain phase has been theoretically proposed HE 
non-invasive techniques have shown that pinning can be 
partially decreased or even increased by applying dc cur¬ 
rents (16( m or ac magnetic fields m- The resulting 
intermediate responses are highly reproducible and inde¬ 
pendent of the previous history. Are these in-between 
responses originated from bulk VL configurations with 
an intermediate degree of disorder? 

To address this question, we performed an experiment 
in a clean NbSe 2 single crystal, combining small angle 
neutron scattering (SANS) with in-situ linear ac suscep¬ 
tibility measurements. In SANS imaging of a perfect tri¬ 
angular VL, first order Bragg peaks (BP) are expected 
to appear as six symmetric sharp spots, matching the 
reciprocal lattice vectors. This still holds for an ideal 
BG, where elastic disorder leads to slow decay of the 
correlations [1], although recent simulations [T8j have 
shown that elastic BG prefracturing would result in a 
slight broadening of the BPs. Conversely, when disor¬ 
der is dominated by dislocations, spatial correlations de¬ 
cay exponentially beyond a characteristic length £ asso¬ 
ciated with the scale at which dislocations appear T9] . 
In such cases BPs have a finite intrinsic width propor¬ 
tional to £ _1 . Experimentally observed BPs are further 
smeared by a convolution kernel due to finite resolution 
[20l . Hence, the degree of disorder can be determined 
from the spread of the scattered intensity, as long as 
£ 1$ L res , a resolution-dependent bound. In a 3D vor¬ 
tex system, the anisotropy of the elastic constants gives 
rise to very different characteristic lengths along the di¬ 
rection of the flux lines (C,l) and in the transverse plane 
(Cl)- Because of the poorer in-plane resolution, in this 
experiment only the longitudinal £/, was determined by 
measuring the integrated intensity on the plane of a single 
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FIG. 1. Scattering geometry: Both the dc field H and the ac 
field h ac are aligned with the c axis of the single crystal, in the 
xz plane. As the whole system is rocked around the y axis, 
the distribution of the scattered intensity I is measured by 
the position sensitive detector, which spans the xy plane, for 
each angle 9 between the field direction and incoming neutron 
beam wave vector K. The maximum intensity is expected at 
Q 0 = Q 0 /K, for the corresponding reciprocal lattice vector 
Qo- 


BP as the sample was rocked through the Bragg condi¬ 
tion (rocking curve/RC). 

The experiment was conducted at the SANS-II beam 
line of the Paul Scherrer Institute’s SINQ facility, using 
the MA11 cryomagnet. The geometry of the experiment 
is sketched in Figure [I] All SANS measurements were 
performed at Tq = 1.9 K. Neutron wavelength was cho¬ 
sen to be An ~ 9 A (±10%) and the applied dc field 
was H = 5 kOe. The incident beam was collimated, 
resulting in a beam divergence ~ ±0.05 deg. The an¬ 
gular resolution along the rocking angle was estimated 
to be w res = (0.107 ± 0.004) deg. A statistical analysis 
showed that curves arising from C,l < L res ~ 40 pm were 
distinguishable from the resolution kernel. RCs wider 
than our resolution were fitted by the convolution of the 
Gaussian resolution with a Lorentzian distribution cor¬ 
responding to an intrinsic A = 2/Q, (see Appendix). 
A large (5.6 x 4.7 x 0.2 mm 3 ) clean NbSe 2 single crystal 
was used for the experiment. The sample is part of a 
batch grown in Bell Labs as described in Ref. EH, with 
T c = (7.115 ± 0.029) K (10 - 90% of the linear ac sus¬ 
ceptibility transition at H = 0). The phase diagram for 
various crystals has been built using a 7-T MPMS XL 
(Quantum Design) and compared with those published 
in the literature P3, showing non significant differences 
(see Appendix). 

Since the goal was to explore the connection between 
structure and dynamics, a non-invasive measurement of 
pinning was required. This was achieved by installing 
a sample holder with mutual inductance coils inside the 
SANS cryomagnet. By adding a small ac magnetic field 
h ac = 2.5 mOe (< 10 -6 IT,/ = 65 kHz), vortices 
were forced to perform small oscillations inside their pin¬ 
ning wells, without modifying their spatial configuration. 
These oscillations propagate through the sample due to 
the vortex-vortex repulsion, with a characteristic pen¬ 
etration depth A qc that is related, in a mean field ap¬ 


proximation [22] . to the effective pinning potential [25] . 
The linearity of the response, as well as the low dissipa¬ 
tion level characteristic of the linear Campbell response 
[251 had been previously verified. Then \ ac and thus the 
effective pinning can be assessed through the in-phase 
component of the linear ac susceptibility x' by means 
of a mutual inductance technique. Lower values of x' 
are associated with stronger pinning, saturating at per¬ 
fect ac shielding (normalized to x! = — 1/47r), whereas x' 
vanishes under complete ac penetration when X ac > sam¬ 
ple size. To probe oscillatory dynamic effects, a shaking 
procedure was applied at certain temperatures: Here 
measurements were paused and a larger sinusoidal field 
with amplitude h s h = 7 Oe, and f s h = 1 kHz was ap¬ 
plied for 1000 cycles before returning to the smaller h ac 
to resume x! measurements. The shaking field amplitude 
was sufficient to induce vortex displacements larger than 
cio throughout the sample. The complete penetration of 
the shaking field into the sample had been previously 
confirmed by measuring the non-linear ac response and 
checking consistency with the irreversible dc magnetiza¬ 
tion determined from M(H) loops (see Appendix). 

Figure [2j a) shows a comparison between y'(T) mea¬ 
sured during a FC procedure and after shaking the VL 
at two different temperatures ± sh below the onset of 
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FIG. 2. (a) Linear ac susceptibility x'(T) measured down 

to To = 1-9 K during a FC cooling procedure (black squares) 
and after shaking the VL at two different temperatures T s ^’ 2 ^ 
(red circles and violet triangles) below the onset of the PE 
( T po ). Shaking temperatures are indicated by vertical lines 
and arrows indicate temperature evolution. Reversibility is 
observed in the warming from To of the VL shaken at 
(dark yellow triangles). The area enclosed by the rectangular 
box is magnified in Fig. |3](a). (b) RC after FC; the black curve 
is the fit of a Lorentzian peak convoluted by the resolution 
kernel, (c) RCs after shaking the VL at (red circles) 
and T s (2) (violet triangles); the curves are Gaussian fits with 
FWHM equal to the instrumental resolution w re s = 0.107 
deg. 
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the PE at T po ~ 6.27 K. There is an unambiguous rise 
in \ ac after shaking, representing a reduction of the ef¬ 
fective pinning. The agreement between both cooling 
x' (T) curves after shaking is consistent with an equili¬ 
bration procedure below T po . Reversibility is observed 
in the warming curve, with y'(T) increasing with T up 
to T po , after which it drops due to the spontaneous pin¬ 
ning enhancement at the PE. In contrast, the absence of 
PE in the linear response during the cooling El is likely 
due to the VL remaining trapped in a metastable highly 
pinned state. This explanation is consistent with the FC 
RC (Figure [2^b)), which is clearly wider than the reso¬ 
lution kernel, in agreement with other experiments |24| . 
Conversely, the two RCs recorded after shaking the VL 
(Figure [2](c)) show resolution-limited BPs (£l > L res ). 
The fact that BPs are well defined (unlike the annulus 
observed for amorphous VLs) indicates that orientational 
order is preserved in all the cases. The broadening of the 
FC configuration RC suggests a dislocated VL with bro¬ 
ken positional order beyond a characteristic volume with 
longitudinal dimension Cl < L res . There is a clear corre¬ 
lation between the ordering of the VL and the decrease 
of the effective pinning. 

While shaking the VL at any T s h < T po leads to an or¬ 
dered configuration, the consequences are strikingly dif¬ 
ferent when the shaking field is applied in the narrow 
transitional region above T po . Figure [3]( a) shows the evo¬ 
lution of y'(T) through a series of processes near the 
transitional region (hatched area) and Figure |4](a) dis¬ 
plays the corresponding RCs. A strongly pinned FC VL 
was cooled to T sh = 6.29 K > T po and shaken, partially 
reducing the effective pinning. The VL was cooled to T 0 
and a narrow RC, still wider than w res was measured, 
indicating a bulk partial ordering. In contrast, applying 
the same shaking burst at the same T ^ (±3 rnK) to an 
initially ordered VL had the opposite effect: The shaking 
field partially increased the effective pinning to a similar 
level and disordered the VL. The RC was broader than 
w res , but narrower than for the maximally disordered FC 
lattice. Because the response after shaking is strongly de¬ 
pendent on the temperature when T s h > T po , the process 
was carefully repeated within A T s h ~ 3 mK, obtaining 
the same qualitative results. To rule out that the par¬ 
tial disordering was of thermal origin, a warming-cooling 
cycle up to T^ h but without applying the shaking burst 
was carried out for a replicated ordered VL. Although the 
cooling x'iT) branch indicates higher pinning than the 
warming branch, it is lower than for VLs shaken above 
T po . Moreover, if there was any associated disorder, it 
was minor and below our SANS resolution. Therefore, 

the “intermediate” degree of disorder found before is only 

(3) 

attributable to the shaking field applied at TV. 

Although the above results are qualitatively evident by 
simply comparing the corresponding intensity distribu¬ 
tion on the position sensitive detector at 6 = do (Figure 



FIG. 3. x'(T) for different thermal and dynamic histories 
(labeled (b)-(e) in panel (a)) and the corresponding detector 
intensity distributions for one VL BP at the Bragg condition 
(panels (b)-(e)). Arrows in (a) indicate temperature evolution 
and the hatched area shows the transitional region. Data from 
Figure [2ja) (black and red lines) is included for comparison. 
The application of a shaking field at T s ® ~ 6.29 K > T po 
results in a unique response (full blue and green symbols) in¬ 
dependent of whether the initial state was a disordered FC VL 
(open blue symbols) or an ordered VL (open green symbols). 
Both intensity distributions (c) and (d) indicate partial disor¬ 
der as compared to the FC intensity distribution (b). Simply 
warming an ordered VL (magenta squares) up to T ~ 
shows moderate pinning enhancement during cooling, but the 
corresponding BP (e) is still sharp and intense. The red boxes 
in (b)-(e) enclose the detector area included in the RCs in Fig¬ 
ure |4] 


UN -(e)), we evaluated the statistical significance of the 
observed differences (see Appendix). Figure |4])b) shows 
the profile likelihood ratio resulting from the fit of 

each experimental RC to one expected from a configura¬ 
tion with an intrinsic width given by £l convoluted by the 
experimental resolution. Correlation lengths obtained af- 

/ o\ 

ter shaking the VL at T^ h ' were L res > £l > £l(FC) at 
the 75% confidence level. Thus, we state that shaking 
the VL in the transitional region results in bulk intrinsic 
configurations with intermediate degree of disorder. 

The total integrated intensity I o (the area under 
the RC) for the different configurations should remain 
roughly constant. Even so, the long-tail decay of the 
structure factor in the BG phase may lead to an under¬ 
estimation of /o [1]. Indeed, for all FC configurations and 
for the configurations obtained after shaking the VL in 
the transitional region we observe a similar Iq. However, 
for the lattice shaken in the BG the integrated intensity 
is slightly, but statistically significant, smaller (see Ap¬ 
pendix). These features support a qualitative difference 
in the positional order decay in the elastic BG configura- 
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FIG. 4. (a) Rocking curves measured after different field- 

temperature histories: FC (black, same as in Fig. §b)), 
after shaking at T® a FC (blue) and an ordered (green) 
VL, and after warming an ordered VL to T = T^ h with¬ 
out shaking at this temperature (magenta). Curves are fits 
of Lorentzian peaks convoluted by the resolution, (b) Cal¬ 
culated profile likelihood ratio Ji? r as a function of Cl up to 
the sample thickness d for each experimental RC. Symbols 
above the curves indicate the maximum („S? r = 1) and er¬ 
ror bars, the 68% confidence interval. For the ordered VL, 
Cl > L res ~ 40 pm ~ 570 a o and is beyond our resolution. 
The profile likelihood of replicate RCs of the FC (black dashed 
line) and the one shaken at T® ± 3 mK (green dashed line) 
are reproducible within the error bars. 


tions and in those obtained after shaking the VL in the 
transitional region, where disorder is probably mediated 
by dislocations. 

In the elastic BG, the in-plane correlation length £j_ 
can be derived from Cl through the tilt and shear elastic 
constants C44 and cqq o 111 EH [26]. From L res , a lower 
limit of 35 a 0 can be estimated for the in plane size of 
the dislocation-free regions corresponding to the ordered 
configurations, which is consistent with results obtained 
in other works in the BG phase [U|23J[27]. How £j_ relates 
to Cl is not as straightforward for metastable VL config¬ 
urations with dislocations. In theoretical works [32. 33]. 
Cl and (j_ are obtained for equilibrium VL configura¬ 
tions from the densities of screw and edge dislocations 
that minimize the free energy; far enough from the amor¬ 


phous limit (i.e., £_l ~ ao ) dislocation densities adapt 
in such a way as to satisfy the elastic ratio. Moreover, 
numerical simulations [ 34 ] and experimental results | 35 ] 
show that Cl monotonically increases with the in-plane 
ordering from the plastic to the elastic regime. In our 
case, the partially disordered VLs are metastable config¬ 
urations (at To) which originated from a dynamic reor¬ 
ganization in the vicinity of the ODT. In this region the 
elastic constants are dispersive and may also be affected 
by fluctuations. Therefore, a precise quantitative esti¬ 
mation of £_l from Cl is not available, but a monotonic 
correlation is still expected (see Appendix for a quanti¬ 
tative estimate). Such a correlation would be compatible 
with an intermediate density of dislocations in partially 
disordered VLs. 

Our results show a clear connection between the linear 
ac response, related to the effective pinning, and the bulk 
spatial correlation of the VL. Configurations having cor¬ 
relation lengths much larger than the lattice parameter, 
but small enough to affect the global response, were ob¬ 
served after shaking the VL in the proposed transitional 
region. These observations support the existence of such 
region and represent the first evidence of a dynamic re¬ 
ordering, driven by oscillatory forces, which results in 
robust bulk configurations with intermediate degrees of 
disorder. A similar behavior was obtained in numerical 
simulations Il2l . where configurations with intermediate 
dislocation densities are accessible from stationary fluc¬ 
tuating dynamic states, as those proposed for colloids 
urn eu- The interplay between these spatial configura¬ 
tions and their associated non-linear dynamics remains a 
fertile ground for further research. 
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APPENDIX 

SANS DATA ANALYSIS 
Resolution model 

We modeled rocking curves (RCs) as the convolution 
of an intrinsic intensity distribution with a smearing res- 
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olution function. The former corresponds to the actual 
g-dependence of the structure factor, while the latter ac¬ 
counts for imperfect collimation and wavelength spread. 
For simplicity, the resolution function was assumed to be 
a Gaussian kernel with standard deviation w res /2. Since 
rocking curves are already integrated within the detector 
plane, only the longitudinal dependence of the structure 
factor is relevant. 

When spatial correlations of the VL decay strongly 
with certain characteristic length scale, the envelope of 
the correlation function can be modeled with an expo¬ 
nential function. So, the structure factor S(qL) (the 
Fourier transform of the correlation function) must be 
Lorentzian-shaped, i.e. 


S(q L ) oc 


1 

1 + (Cl9l) 2 ’ 


(1) 


where Ql is the corresponding correlation length. There¬ 
fore, we took the intrinsic distribution to be the Cauchy- 
Lorentz distribution centred at the Bragg angle 9q , mul¬ 
tiplied by a proportionality constant Iq. Then, its half 
width at half maximum (HWHM) 7 is directly related to 
Cl by 


1 

Q^C~l 


( 2 ) 


where Qo is the magnitude of one reciprocal vector of the 
VL. The resulting intensity profile is 


1(d) = 


7 


2 e 


- 2 (^) S 


'-oo 7 T 7 2 + (O' - do) 2 V'StT 


7 TW r 


-dfl'. (3) 


In the limiting cases of long correlation lengths QoCl 
l/w res , (i.e. when C,l is beyond our resolution) the in¬ 
trinsic distribution can be replaced by a Dirac delta, so 


m 


2 /n 


\f2/Kh 


■ exp 


-2 


'-0c 


(4) 


and the experimental RCs can be just modeled by the 
Gaussian resolution distribution. 


Data reduction 

Neutron count data were preprocessed by the software 
GRASP [35]. Detector counts were normalized by mon¬ 
itor counts. Then, after background subtraction, inten¬ 
sity integrated over a fixed sector of the detector was 
calculated for each image, and resulting rocking curves 
(■ 0j,lj ), j = were exported, together with their 

error estimates A Ij. 


Data analysis 

Statistical analysis was carried out using the software 
Octave E3- It consisted of three stages: 


1 . preliminary analysis, 

2 . resolution estimation, 

3. intrinsic widths estimation. 


Preliminary analysis 

Non linear optimization often requires set of initial pa¬ 
rameters close enough to the optimal in order to work 
properly. The aim of this stage was to obtain suitable 
estimates of the parameters for each rocking curve indi¬ 
vidually, to be used in later stages. We did this by fitting 
our convolution model individually to each rocking curve. 

Because of the convolution, the parameters w res and 
7 are expected to be highly correlated, making it diffi¬ 
cult to fit both simultaneously. To overcome this issue, a 
rough prior theoretical estimate of w res , based on the ex¬ 
perimental geometry [33], was calculated. This resulted 
in w res = (0.104 ± 0.003) deg. A modified likelihood 
function Jz? was then constructed and parameters were 
obtained from minimizing 


-21ogjSf(/ 0 , 0 O ,7, Wres) = 

2 N 

+E 


Wres 


Aw, 


3 =1 


(ij-m) 

\ A/, 


2 


(5) 


Resolution estimation 

The preliminary analysis showed that configurations 
shaken at low temperatures (examples are shown with red 
and violet symbols in Figure[2](c)) were resolution limited 
(7 -C w res ). Therefore, we used the long-correlation- 
length approximation (Eq. |4j) to model these rocking 
curves. In turn, it enabled us to make an independent 
estimation of the resolution based solely on neutron data. 

The resolution profile likelihood Jz? ' W (w res ) was con¬ 
structed by maximizing the global likelihood function 
(i. e., including all rocking curves), constrained to w res 
taking a series of fixed values. We obtained w res = 
(0.107 ± 0.004) deg, which is consistent with our prior 
estimate. This also revealed that the dependence of the 
estimated 7 on the resolution width was common to all 
the curves. Examples of the resulting fitting curves are 
shown in full red and violet lines in Figure [2jc). 

Next, we determined the minimum intrinsic width 7 
which could be eventually distinguishable from the zero- 
width approximation. This was done by generating ran¬ 
domized rocking curves of different widths with the con¬ 
volution model, using the above estimation of w res . Sim¬ 
ilar number of points and relative errors as in our ex¬ 
perimental data were considered. After fitting Gaussian 
functions to generated data, a threshold of > 90% out 
of 200 replicates resulting in non-overlapping error bars 
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(one standard deviation) between w res and the fitted 
width was applied. This resulted in 7 > 0.014 deg, or 
Cl < 40 pm ~ L res (vertical dashed line in Figure |4](b)). 

Intrinsic correlation length estimation 

Lastly, the profile likelihood of 7 for each rocking curve 
was constructed by fitting only Iq and 6q, while letting 
w res be the estimated above (fixed). This procedure re¬ 
duced marginal errors in 7 by ignoring the covariance 
with w res . Since a change in w res would only shift all 
estimated 7 together, this approach is better suited for 
comparing different configurations instead of finding the 
exact correlation length. The resulting likelihood ra¬ 
tios Jzf r (Ci) (curves in Figure |4|b)) where then obtained 
through the relationship between and 7 (Eq. [2|. 

As long as the elastic approximation is valid, the in¬ 
plane characteristic lengths £_l in equilibrium configura¬ 
tions can be estimated from the elastic ratio 

Cl ~ Cl \jCQs/ca, (6) 

where C66 and C44 are the shear and tilt elastic VL moduli 
respectively. In the local approximation [26] , 

C 44 /C 66 ~ 16tt (A/a 0 ) 2 , (7) 

where A is the London penetration depth and do is the 
mean intervortex distance. In our case 1251 . A (T = 2 
K)/ao ~ 2 , so the corresponding in-plane characteris¬ 
tic lengths in the BG phase are expected to be C±bg > 
L res /14: ~ 35 do- On the other hand, the fact that 
all the BPs are resolution limited, indicates that all the 
measured VL configurations are not in the amorphous 
limit. 

Far from the amorphous limit, the elastic approxima¬ 
tion is still roughly valid for equilibrium configurations 
with dislocations [32112$ • However, as mentioned in the 
main text, the partially disordered VLs are metastable 
configurations at Tq = 2A', originated from a dynamic 
reorganization in the vicinity of the ODT, and there¬ 
fore a quantitative calculation of C_ from (L is not avail¬ 
able. A rough estimation using the elastic limit at T 0 
(where the local approximation for the elastic constant 
is valid) gives Cl(FC) ~ 5 do and Cl ~ 12 do for con¬ 
figurations with intermediate disorder. A better esti¬ 
mate can be made in the framework of a defective VL 
[32] I33| at high temperature (where the dynamic reor¬ 
ganization occurs). In this model each correlated vol¬ 
ume ClCI is limited by dislocations spaced in Ci_ and 
in the transverse and longitudinal directions respec¬ 
tively. To calculate Cl (Cl) we follow Eqs. (9) and ( 10 ) 
of Ref. (321 (where the notation is D c = L c = Cl) 
and a non local approximation for the tilt modulus [25] 
c 44 (Cl) ~ c 44 (0)/(A 2 fc 2 + 1) ~ C4 4 (0)/(27tA/C_l) 2 where 
the penetration depth A £_l/27t. The resulting Cl (Cl) 


are A-independent in this limit and the in-plane correla¬ 
tion lengths corresponding to the experimental Cl result 
C_l(FC) ~ 9 do and Cl ~ 18 do for intermediate disor¬ 
dered configurations, slightly larger to that obtained in 
the elastic limits. Although these values are not quanti¬ 
tatively exact, these estimates could be useful to design 
future SANS experiments. 

A comment on total integrated intensities 

Equation ([ 3 ]) is properly normalized so that the pa¬ 
rameter Iq represents the total integrated intensity. The 
maximum-likelihood estimate for I 0 and the correspond¬ 
ing variance were obtained along with the other param¬ 
eters discussed above for each RC. In order to compare 
integrated intensities corresponding to different config¬ 
urations, the dispersion in the intensities fitted to RCs 
measured after different histories were statistically as¬ 
sessed by means of \ 2 difference tests of nested mod¬ 
els. The results of the tests are presented indicating the 
statistic x 2 , the degrees of freedom df , and the corre¬ 
sponding p- value (i. e., the probability of getting a x 2 
greater than the observed one solely due to randomness, 
under the assumption that a common mean is shared). 

First we note that the set of all the values of 7 0 had far 
more dispersion (relative to their errors) than expected 
for a single mean value model (x 2 = 42.48, df = 10, 
p < 0.001). In order to assess whether this excessive 
dispersion was connected to the different histories, the 
dispersion within sets of Iq corresponding to each pro¬ 
cess was evaluated. All intensities corresponding to FC 
configurations were indistinguishable (x 2 = 1.37, df = 2, 
p = 0.505) and intensities corresponding to intermediate- 
order configurations were only marginally distinguish¬ 
able (x 2 = 6.15, df = 2, p = 0.046). Moreover, merg¬ 
ing the sets of intensities for FC and intermediate-order 
configurations did not increase the dispersion enough to 
reject a single mean intensity for both (Ay 2 < 0 . 01 , 
A df = 1, p = 0.938). On the contrary, BG configu¬ 
rations presented a higher degree of variability, beyond 
their estimated variance (x 2 = 14.78, df = 4, p = 0.005) 
and the set of intensities for BG configurations could 
not be merged with FC and intermediate-order config¬ 
urations without increasing the dispersion significantly 

(Ax 2 = 20.18, A df = 1 ,p< 0.001). 

These results suggest that all RCs for FC and 
intermediate-order configurations share essentially the 
same value for the parameter I 0 . Instead, fits to RCs 
for BG configurations yielded statistically distinguishable 
values for the total integrated intensity, which averaged 
about 19% below the mean I 0 for FC and intermediate- 
order configurations. As mentioned in the main text, an 
underestimation of I 0 is likely to result from the long-tail 
decay of the structure factor in the BG phase, further 
supporting the observed structural differences. 
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VORTEX PHASE DIAGRAM 

In this supplementary section, the vortex phase dia¬ 
gram corresponding to our samples is presented and some 
features are compared with those reported in the litera¬ 
ture. The magnetic field and temperature ranges used in 
our experiments are also included. As the sample used in 
the SANS experiment (in the following, sample A) was 
too large to be measured in the 7-T MPMS XL Quan¬ 
tum Design or in our cryostats, another crystal from the 
same source (sample B, V ~ 0.25 x 0.25 x 0.02 cm 3 ) and 
a smaller sample obtained by cutting an edge of sample 
A (sample C, V ~ 0.1 x 0.1 x 0.02 cnr 3 ) were charac¬ 
terized instead. Beyond the expected differences due to 
the larger area of sample B, results obtained from both 
samples are identical. In Figure [5] the phase diagram 
obtained from non-linear ac susceptibility curves x'(H) 
(open symbols) and x'(T) (full symbols) measured in 
sample B is shown in panel (a), including the transitional 
region (hatched area in Fig. [5]). The dark yellow line is 
the spinodal, lying above the transitional region, that de¬ 
marcates the limit of stability with respect to fluctuations 
toward the thermodynamically stable state, extracted 
from Ref. [T7| and rescaled using the critical parameters 
of our sample T c = 7.11 K and H c 2 ( 0 ) = 54.45 kOe (esti¬ 
mated by extrapolating the line x'(H- T) = 0 to T = 0). 
This line matches the abrupt change in slope observed in 
both x'(T) ( T sp in panel (b)) and x'(-ff) ( H sp in panel 



FIG. 5. (a) Vortex phase diagram obtained from x'(T) (full 
symbols, example in panel (b)) and x'(H) (open symbols, 
example in panel (c)). Below H ~ H po (T) a shaking field 
orders the VL. The hatched region indicates the transitional 
region, where we claim that the VL dynamically reorganizes 
in configurations with intermediate degree of disorder. Above 
H > H P (T) there is a narrow region where disordered VLs 
are stable. No history effects are observed beyond H sp (T), 
which matches the spinodal line described in Ref. [17] (dark 
yellow line). 


(c)), above which no history effects are observed. The 
fact that our experiments are performed below the spin¬ 
odal line indicates that fluctuations are not dominant in 
the underlying physics. However, due to the close vicin¬ 
ity, some small influence cannot be discarded. The dotted 
horizontal line indicates the magnetic field applied in our 
SANS experiment. 

MAGNETIC FIELD PENETRATION AND 
VORTEX DISPLACEMENTS 

The full bulk penetration of the shaking fields, with the 
consequent vortex displacements in distances larger than 
ao, is essential to justify our claim of a bulk dynamic 
reordering. The modification of the bulk VL structure 
factor observed in the SANS experiment is the main di¬ 
rect evidence of such a reordering. In this supplementary 
section, we present some numerical estimations, together 
with additional experimental results obtained by magne¬ 
tization and non-linear ac susceptibility measurements, 
that provide a consistent scenario regarding this crucial 
issue. 

Estimation of vortex displacements 

A simple calculation provides a rough estimation of 
the mean vortex displacement associated with the shak¬ 
ing field penetration. Under full ac field penetration, the 
maximum variation of the VL parameter a 0 during a cy¬ 
cle of the shaking field is 

Aao ~ {a 0 /2)(h sh /B) ~ (a 0 /2 )(h sh /H). 

In our SANS experiments, h s h ~ 7 Oe and H = 5 kOe, so 
Aao ~ 0.05 nrn. Therefore, near the centre of the sample, 
vortices only need to move short distances u > Ado to 
accommodate to the change in the density. However, the 
propagation of the density change throughout the sam¬ 
ple implies that the mean displacement increases with 
the distance r from the centre as u(r) ~ r (Ado/do) . 
Consequently, in our experiments we expect u(r ) > do 
for r > 100 pm, reaching u[r) as large as several microns 
near the edge of the sample (r ~ 2.5 mm). In conclu¬ 
sion, during the shaking procedure, most vortices move 
long distances u > ao, allowing a dynamic reorganiza¬ 
tion. On the other hand, vortices that enter and exit the 
sample crossing its edge during the shaking procedure 
are present only in a thin micrometric annulus near the 
sample boundary (less than 2% of the sample area) and 
may only marginally affect the average bulk structure. 

The claim that the small ac field h ac = 2.5 mOe, ap¬ 
plied during the linear ac susceptibility measurements, 
does not modify the VL configuration, i.e. the technique 
is non-invasive, has been extensively discussed in previ¬ 
ous publications m-- the linear non-dissipative Camp- 










bell regime j23] holds when each vortex performs small 
harmonic oscillations around a minimum of the effective 
pinning potential well. These small oscillations propa¬ 
gate inside the sample up to the Campbell penetration 
depth, generally smaller than the sample dimension. The 
same reasons outlined above (replacing h s h ~ 7 Oe by 
h ac = 2.5 mOe) lead to the conclusion that during the ac 
measurement even outmost vortices move typically less 
than 3 nm < £ 0 < a 0 ; where £o is the coherence length, 
not modifying the VL configuration. 


Linear and non-linear ac regimes: ac field 
penetration and dissipation 

The full penetration of the shaking field was previ¬ 
ously confirmed by measuring the non-linear ac response 
in sample C, and results are shown is Figure [d] (color 
curves): Below T pi the normalized ac susceptibility is am¬ 
plitude dependent, with a maximum out of phase compo¬ 
nent Xma. x ~ 0.24, characteristic of the ac Bean Critical 
response .32]. At T > 5 K, an ac field with amplitude 
h ac = 3 Oe (the maximum available in the MPMS) and 
frequency / = 1 kHz fully penetrates the sample (%' ~ 0); 
an even higher penetration is expected for the shaking 
field h s h ~ 7 Oe. Notice that, despite the smaller volume 
of sample C, because samples A and C have exactly the 
same thickness d, similar non linear ac field penetration 
will hold [391. In the same figure, the linear (i.e. ampli¬ 
tude independent) ac susceptibility, measured in sample 
A with the SANS setup, is also included. Below T p , the 
dissipative out-of-phase component ~ 0, as expected 
in the linear Campbell regime mm- 


Magnetic field penetration in dc measurements 

The conclusions of the above paragraph are also con¬ 
sistent with results obtained by magnetization measure¬ 
ments. Figure [7] shows examples of magnetization curves 
at low fields H measured at different temperatures. The 
shape of the M(H) loops is characteristic of low pin¬ 
ning thin superconductors in a transverse geometry m, 
where geometrical barriers inhibit the magnetic field en¬ 
trance up to a field H en higher than the effective lower 
critical field H*^ (~ 0.15 H c \ in our experimental geom¬ 
etry). At fields H H en , geometrical barriers do not 
play any significant role but a small irreversible magneti¬ 
zation M rrr (H). arising from weak bulk pinning, persists. 
It can be observed in the inset of figure [7] that M irr de¬ 
cays quickly with increasing field. At H = 5 kOe (the 
magnetic field applied in our SANS experiments), the low 
irreversibility (47r Mi rr < 1 Oe) is consistent with the full 
bulk penetration of the shaking field h s h ~ 7 Oe. 
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FIG. 6. Linear (black symbols) and non-linear (color symbols) 
ac response for different ac amplitudes, measured at H = 5 
kOe. The sharp zero-field transition is also included in light 
grey symbols. In the linear response ( h ac = 2.5 mOe) x" ~ 0 
below T p , consistently with a non-invasive Campbell regime. 
The highest-amplitude ac field fully penetrates the sample, 
for T > 5 K. Arrows in the linear curves show the direction 
of temperature variations. Non linear curves were recorded in 
cooling procedures. 



FIG. 7. (a) Magnetization loops M(H) at low fields, for sev¬ 
eral temperatures. The shape is characteristic of low pinning 
thin superconductors in a transverse geometry; geometrical 
barriers inhibit the magnetic field entry up to H en > H* ± 
(the effective lower critical field). Inset: At H H en , a 
small irreversible magnetization Mi rr (H), arising from weak 
bulk pinning, persists. At H = 5 k Oe, 4irMi rr < 1 Oe, 
consistently with the full bulk penetration of h s h ~ 7 Oe. 
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